/*=========================================================================
 Copyright (c) 2008, Nokia Corporation and/or its subsidiary(-ies).
 All rights reserved.

 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
 
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Nokia Corporation and/or its subsidiary(-ies)nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.

 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
 MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
 THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
 OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR 
 TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=========================================================================================*/

#ifndef SVD_H
#define SVD_H
#include <e32base.h>
#include "ncvMatrix.h"
#include "ncvVector.h"
#include "ncvFixed.h"

/**
 * \class CSVD
 * \brief "Economy size" singular value decomposition of a given real matrix.
 *
 *  <pre>
 *  Usage:
 *  CNokiaCVMatrix *a = CNokiaCVMatrix::NewL(rows,cols);
 *  CSVD *svd;
 *  TRAPD(err, svd = CSVD::NewL(*a))
 *  ... //Handle possible errors
 *  CNokiaCVMatrix *v = svd->V();
 *  ...
 *  delete a;
 *  delete svd;
 *  delete v;
 *  </pre>
 *  
 * The algorithm produces the "economy size" decomposition. 
 * If  X  is m-by-n with m > n, then only the first n 
 * columns of  U  are computed and  S  is n-by-n.
 *
 * Householder bidiagonalization and a variant of the QR algorithm are used
 *
 * The SVD algorithm is derived from JAMA/TNT, the Java Matrix Package.
 * /see http://math.nist.gov/javanumerics/jama/ for more information
 * More information about the implementation can be found e.g. from
 * George E. Forsythe et al. from Computer Methods for Mathematical Computations.
 */

/*
 * The max number of iterations to be made during the main iteration of the svd
 * If exceeded the system will throw KNokiaCVErrToManyIterations
 */
#define MAX_ITERATIONS 30

class CSVD: public CBase
	{
private:
	CSVD();
	/*!
	 * \brief Constructs the SVD
	 * \param aMatrix the matrix to decompose
	 * \param aGenU Tells whether to generate the U matrix or not (ETrue => Generated)
	 * \param aGenV Tells whether to generate the V matrix or not (ETrue => Generated)
	 */
	void ConstructL(CNokiaCVMatrix &aMatrix, TBool aGenU = ETrue, TBool aGenV = ETrue);

	/*!
	 * \brief Does the singular value decomposition A = U*S*V'of the
	 *        matrix given in the constructor. U and V are generated by default, but
	 *		  not if aGenU and aGenV are set to EFalse
	 *
	 *  \note KNokiaCVErrToManyIterations from ncvErrors.h will be thrown if the
	 *  algorithim exceeds the defined MAX_ITERATIONS iterations 
	 *  and KErrNoMemory if we run out of memory
	 *
	 *  \note the function can leave
	 */
	void Decompose();
	
	/*!
	 * \brief Bidiagonalize the matrix given in the constructor
	 *
	 * Forms the Householder bidiagonalization A = XBY^t 
	 * of the matrix given in the constructor.
	 * 
	 *					  |s e  |
	 * where e.g. 3x3 B = |  s e|
	 *			          |    s|
	 *
     * The diagonal elements are stored in the array iS and the 
	 * super-diagonal elements are stored in the array iE.
	 *
	 * \note Householder bidiagonalization is not so good with sparse matrices
	 *		 but excellent for dense matrices.
	 */
	void Bidiagonalize();
	
	/*!
	 * \brief Generates the U matrix of the bidiagonalized matrix
	 */
	void GenerateU();
	
	/*!
	 * \brief Generates the V of the bidiagonalized matrix
	 */
	void GenerateV();

	/*!
	 * \brief Helper method calculates hypotenuse between two points
	 * \return The hypotenuse between the given numbers as TReal
	 */
	TReal Pythag(TReal a, TReal b);


public:
	IMPORT_C ~CSVD();

	/*!
	 * \brief Creates the singular value decomposition of the given matrix
	 * \param aMatrix the matrix to decompose
	 * \param aGenU Tells whether to generate the U matrix or not (ETrue => Generated)
	 * \param aGenV Tells whether to generate the V matrix or not (ETrue => Generated)
	 *  KNokiaCVErrToManyIterations (ncvErrors.h) will be thrown if the
	 *  algorithim exceed 30 iterations and KErrNoMemory if we run out of memory
	 * 
	 * \note the function can leave
	 * 
	 */
	IMPORT_C static CSVD* NewL(CNokiaCVMatrix &aMatrix, TBool aGenU = ETrue, TBool aGenV = ETrue);

	/*!
	 * \brief Getter for the V generated by the SVD
	 * \return the generated V
	 */
	IMPORT_C CNokiaCVMatrix* V();
	/*!
	 * \brief Getter for the V' generated by the SVD
	 * \return the generated V transpose
	 * \note the function can leave
	 */
	IMPORT_C CNokiaCVMatrix* TransposeV();
	/*!
	 * \brief Getter for the U generated by the SVD
	 * \return the generated U
	 */
	IMPORT_C CNokiaCVMatrix* U();
	/*!
	 * \brief Getter for the calculated singular values by the SVD
	 * \return the singular values in a diagonal matrix
	 */
	IMPORT_C CNokiaCVMatrix* S();
	/*!
	 * \brief Getter for the calculated singular values by the SVD
	 * \note the function can leave
	 * \return the singular values in a CFixed array in order s[0] >= s[1] >= ... >= s[iN]
	 */
	IMPORT_C CFixed* Singularvalues();
	
	/*!
	 * \brief Return the Norm2 of the given matrix = the largest singular value
	 * \return The largest singular value
	 */
	IMPORT_C CFixed Norm();


private:
	
	/*!
	 * \brief Contains the orthogonal column vectors of the decompositition 
	 * (right side U) iGenU is set to true, otherwise used as a temporary storage
	 */
	TReal** iU;
	
	/*!
	 * \brief Contains the given matrix in TReal format. Not altered during computation.
	 */
	TReal** iA;
	
	
	/*!
	 * \brief Contains the orthogonal vectors of the decompositation 
	 * (left side V) iGenV is set to true, otherwise used as a temporary storage
	 */
	TReal** iV;

	/*!
	 * \brief A temporary storage
	 */
	TReal* iW;

	/*!
	 * \brief The super diagonal elements of the bidiagonalized matrix
	 * are stored in iE.
	 */
	TReal* iE;
	/*!
	 * \brief Contains the actual singular values
	 * 
	 * During computation the diagonal elements of the bidiagonalized matrix
	 * are stored in iS.
	 */
	TReal* iS;

	/*!
	 * \brief Number of columns of the given matrix
	 */ 
	TInt iN;

	/*!
	 * \brief Number of rows of the given matrix
	 */ 
	TInt iM;

	/*!
	 * \brief Min(iM-1, iN), used during computation
	 */ 
	TInt nct;
	
	/*!
	 * \brief Max(0, Min(iN-2,iM)), used during computation
	 */ 
	TInt nrt;

	/*!
	 * \brief Min(iM, iN)
	 */ 
	TInt nu;

	/*!
	 * \brief set to ETrue to compute U
	 */ 
	TBool iGenU;

	/*!
	* \brief set to ETrue to compute V
	 */
	TBool iGenV;
	
	
	};

#endif // SVD_H
